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We use a self-consistent Hartree-Fock approximation with realistic Coulomb interactions for 7T-band electrons 
to explore the possibility of broken symmetry states in weakly disordered ABC stacked trilayer graphene. The 
competition between gapped and gapless broken symmetry states, and normal states is studied by comparing 
total energies. We find that gapped states are favored and that, unlike the bilayer case, gapless nematic broken 
symmetry states are not metastable. Among the gapped states the layer antiferromagnetic state is favored over 
anomalous and spin Hall states. 
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I. INTRODUCTION 

The electronic structure of few layer graphene 1 2 systems 
consists of pairs of bands which cross, or narrowly avoid 
crossing, near the Fermi level. Because the number of band 
pairs depends in an interesting way on the stacking arrange- 
ment, and because both semi-metallic and semi-conducting 
behaviors occur, this family of two-dimensional materials pro- 
vides an attractive playground for the study of electron inter- 
action effects in systems 34 with approximate Fermi points. 
For example, interactions lead to a marginal Fermi liquid 
behavior in neutral single layer graphenepE3 and to broken- 
symmetry ordered states in bilayers.^"^! 

The recent surge of interest in ABC stacked trilayer 
graphene^SK3 motivates theoretical studies of the electron in- 
teraction induced instabilities that are expected when these 
structures are weakly disordered. Instabilities are favored in 
ABC trilayers by extremely flat crossing 2 of a single-pair of 
bands at the neutral system Fermi level, and by exchange en- 
ergy frustration associated with mometum- space textures in 
the valence band wavef unctions. 28 In the trilayer case the 
competition between competing broken symmetry states is 
massaged by weak remote neighbor hopping processes that 
reshape the bands at energies within ~ 20 meV of the cross- 
ing point 29 This energy scale should be compared to the ~ 1 
meV scale of analogous processes in bilayer graphene®!. The 
remote hopping processes are therefore more likely to play 
a prominent role in determining how the system responds to 
electron-electron interactions in the trilayer case. 

The broken symmetry states that have been discussed in the 
bilayer graphene literature can broadly be classified either as 
gapped phases with broken layer inversion symmetryP"^ or 
as gapless nematic states that lower rotational symmetry. 11 Al- 
though the two types of states in principle should have clear 
experimental signatures, it has not yet been possibl e 13 " 15 ! to 
achieve a universally accepted consensus on the character of 
the ground state because of the complicating role of resid- 
ual disorder. Studies of ABC stacked trilayer graphene could 
prove to be more unambiguous because its bands are flatter 
and interaction effects correspondingly stronger, while disor- 
der strengths should be comparable. 

In this paper we present a study of the competition between 
gapful and gapless states in ABC stacked trilayer graphene, 
including the effects of weak remote-hopping processes which 



can dominate band dispersion very close to the band-crossing 
(Dirac) point. Our study is based on a six-band ^-orbital 
tight-binding model, combined with long-range Coulomb in- 
teractions treated using a Hartree-Fock mean-field-theory. 
The quasiparticle band-structures of both gapped and gap- 
less states are reshaped when interactions are included. We 
find that gapped phases are favored over a wide range of the 
hopping-parameter model space. In mean-field theory the en- 
ergy difference between gapped and gapless states is typi- 
cally smaller than ~ 10 -7 eV per carbon atom. The small 
condensation energy reflects the fact that only single-electron 
states close to the band crossing points participate in order- 
ing. The strength of the direct hopping process between low- 
energy sites on the outer layers of the trilayer, 72, plays an 
especially important role in determining the character of the 
ground state. 



II. MODEL HAMILTONIAN 

We describe ABC trilayer graphene using a lattice model 
Hamiltonian with one atomic 2p z orbital per carbon site. We 
label the six sublattice sites illustrated in Fig. 1(a) A, B, A 7 , 
B f , A" , B"\ the A and B" sites avoid near-neighbor inter-layer 
coupling and for this reason are low-energy sites which are 
dominantly occupied by electrons close to the band crossing 
points. With this convention, the six band tight-binding model 
Hamiltonian of ABC trilayer graphene is: 
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with a = 2.46A using the same triangular lattice vector con- 
vention as in Ref . 14181 . The global minus sign in front of the 
Hamiltonian means that ^-bonding bands have lower energy 
than anti-bonding bands when the 7 parameters are positive. 
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FIG. 1 : Trilayer graphene unit cell and 7T-band structure, (a) Schematic representation of the hopping processes included in our model, (b) 
Band structure near the Dirac point K = (4n/3a,0) for different values of the remote hopping parameters. The upper row shows 2D band 
structures seen from a view rotated by 30° with respect to vertical while the lower row shows the same information expressed in terms of 
contour plots. When non-zero, the hopping parameters have the values 72 = 0.01 eV, 73 = 0.3 eV. We have set 74 = 75 = throughout our 
calculations. Wave vectors k are in units of a~ l . The 72 term splits the Brillouin-zone corner cubic band crossing into three Dirac cones 
(linear band crossings) located at the vertices of an equilateral triangle. The trigonal 73 term acting on its own results in four band-crossing 
points, including one at the Brillouin-zone corner. When both terms are present simultaneously the gapless point at K disappears. When the 
signs of both terms are the same they tend to produce opposing triangular distortions, whereas if they have opposite signs their triangular 
distortions are reinforcing. The orientation of the triangular distortion for each of the above parameters is sign dependent. In realistic band 
structures the 72 term dominates over 73 resulting in three Fermi points with linear dispersion. The nematic term Tv captures the influence of 
a layer relative- sliding sliding strain and breaks the triangular rotational symmetry of the bands. We used the value = 0.02 eV in the above 
illustration. 



In most of our calculations we have used graphite hopping 
parameter values which are similar to those in Ref. l3T1l : 
70 = 3.12 eV, 71 = 0.377 eV, 72 = 0.01 eV, 73 = 0.3 eV. We 
specifically address the importance of the signs of the remote 
72 and 73 hopping parameters. The near-neighbor intralayer 
and interlayer hopping processes 70 and J\ are responsible for 
broad features of the band structure, while the 72 and 73 pa- 
rameters have their main impact close to the band-crossing 
points. This model qualitatively reproduces the ab initio band 
structure in Ref. l32l . in particular capturing the orientation 
of the triangle formed by the three band-crossing points close 
to the Brillouin-zone corner. We have ignored the ABC tri- 
layer 74 and 75 processes that break particle-hole symmetry, 
and other small onsite terms that are often introduced in mod- 
els of graphite, because they do not visibly alter the low en- 
ergy features of the bands in ABC trilayer graphene. 

Using a model similar to that used previously for bilayer 
graphene! 34 l 35 i we have also examined the influence of a 
term in the Hamiltonian that is intended to capture the influ- 
ence on low-energy states of an interlayer relative-translation 



/k r ), introducing 



strain. We write 7^ = ^exp(— k — 

a damping factor which makes the term small away from 
the Brillouin-zone corners, where this form for the strain 
Hamiltonian becomes inaccurate, by setting k r = Yi/hVp = 
0.0573A- 1 . 

Because there is some confusion in the literature on the 
signs of the remote hopping processes, we have also consid- 
ered other sign choices for 72 and 73. As shown in Fig. 1(b) 
direct hopping 72 between the low energy sites A, B'" gives 
rise to three Fermi points at the vertices of a triangle centered 
on the Brillouin-zone corner. The trigonal warping (73) pro- 
cess which connect the A, B' and A 7 , B" sites is also respon- 
sible for a trigonal distortion that leads to four Fermi points 



near K, as in bilayer graphene. Each one of the three Fermi 
points contribute to a phase winding of 2tc for a total 6n phase 
winding along paths that encircle all three points, as expected 
in ABC trilayer graphene 2 . (We use the term Fermi point to 
refer to a band crossing that is tied to the Fermi level of a 
neutral ABC trilayer. The band crossing points are exactly at 
the Fermi level because we have neglected particle-hole sym- 
metry breaking terms in the band structure model.) Both 72 
and 73 terms break circular symmetry near the Dirac point 
by splitting a single Fermi point with cubic band dispersion 
into multiple Fermi points with linear dispersion. The orien- 
tations of the triangular distortion due to 72 and 73 are opposite 
when both hopping parameters have the same sign. First prin- 
ciples band structure calculations suggest that 72 dominates 
over 73 and determines the shape of the bands near the Dirac 
point. (Note that 72 has a much greater influence on the two 
low-energy states than 73 for a given numerical value because 
it couples them directly, whereas 73 acts virtually via high- 
energy states.) When both terms are present simultaneously 
and have the same sign the band structure can have a a hy- 
brid shape; for some parameters values the bands consist of 
two intertwined triangles with opposite orientations that can 
exhibit up to nine Dirac cones. The additional parameter, Tv 
couples A, B' and A 7 , B" , like the 73 term, but without an ac- 
companying factor /(k). The Tv term qualitatively describes 
a band deformation that lower the crystal rotational symmetry, 
and is similar to the model used to mimic a small layer- sliding 
structural deformation in bilayer graphene 35 . This term is also 
useful to seed lowered rotational symmetry gapless states in 
our Hartree-Fock calculations. 



Electron-electron interaction effects are treated in an unre- 
stricted Hartree-Fock approximation™ which allows symme- 



tries to be broken: 
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where c^, c k ^ are Bloch state creation and annihilation op- 
erators, and A = (l,o) lumps lattice and spin indices. The 
Hartree and Exchange Coulomb integrals in Eq.([3]>, 
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involve sums over reciprocal lattice vectors G. In these equa- 
tions S/ is the (2D projection of the) position of the sublat- 
tice in the unit cell. We used an isotropic atomic orbital form 



factor /(q) = Jdre^ ft (r)| z = (1 - {r qf)/{{\ + (r q) 2 ) 4 ) 



with an artificially large atomic radius r Q = 3a o /y30 to ac- 
count for sp2 orbital polarization. 4 Here a = a/ (2y/3) is 
the covalent bond radius of carbon. The two-dimensional 
Coulomb interactions in Eqs. are defined by V 11 ' (q) = 
lite 1 1 (|q|e r ) when the sublattice indices / and /' refer to 
the atoms in the same layer, and (lite 2 / (|q|e r )) exp [— |q| d] 
when they refer to atoms in layers separated by a distance d. 

We used an effective dielectric constant £ r = 4 in our calcu- 
lations, partly to account for dielectric screening by surround- 
ing material and partly to account for the well known tendency 
of Hartree-Fock approximation, which neglects screening, to 
overestimate exchange interaction effects. 33 The present im- 
plementation of the lattice model Hartree-Fock mean-field 
theory follows closely the method described in Refs. 14181 
for single and bilayer graphene which also used a momen- 
tum space representation of the Coulomb interaction. One 
difference in the present implementation is that we sample 
the full Brillouin zone without taking advantage of the sym- 
metry of the crystal in order to allow for the possibility of 
broken rotational symmetry nematic phases. Because of the 
greater importance of states near the Dirac point we have sam- 
pled momentum space non-uniformly; for k-points closer than 
~ 0.5 /a to the Dirac point (where a = 2.46A is the triangular 
lattice constant of graphene) we have used a sampling density 
corresponding to 2304 x 2304 points in the entire Brillouin 
zone. Outside this region we used a matched but coarser sam- 
pling with density corresponding to 18 x 18 points in the Bril- 
louin zone. For a given sampling density, the Hartree-Fock 
equations are solved iteratively and converged to ~ 10 -11 eV 
per carbon atom in total energy. 



III. GAPPED AND GAPLESS STATES 

As in the AB bilayer case, the low energy valence band 
states of ABC graphene are given approximately by equal 
weight coherent sums of top and bottom layer wavefunctions 
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FIG. 2: Self-consistent Hartree-Fock calculations iterated from 
seeds for the gapped and gapless nematic states in ABC trilayer 
graphene. The gapped solution (a) has been obtained starting from 
the layer antiferromagnetic initial condition, while the gapless self- 
consistent solution (b) has been obtained by seeding with a nematic 
Yn term. Note that the gapless solutions restore the rotational sym- 
metry of the crystal lattice that was broken by the Yn term, indi- 
cating that nematic order is not stable at the mean-field level in the 
trilayer case. When the remote hopping parameters 72 = 0.01 eV and 
73 =0.3 eV are accounted for they induce a triangular distortion of 
the bands near the Dirac point and determine the angles at which the 
band crossings occur. These processes do not have a large influence 
on the gapped ground state. 



with momentum-dependent phase differences. The gapped 
broken symmetry states spontaneously increase weight in one 
of the two layers, whereat the nematic states break the lat- 
tice rotational symmetry of the inter-layer phases. In the 
following we present the results of our 7T-band Coulomb- 
interaction Hartree-Fock study. This mean-field-theory cal- 
culation we perform cannot be fully quantitative because it 
accounts for screening in an ad-hoc way which might be quan- 
titatively inaccurate, and because it neglects higher-order cor- 
relation effects. We believe though that our results provide 
some insight into the competition between different potential 
ordered states, and in particular into the way this competition 
is influenced by band structure features particular to ABC tri- 
layer graphene. We first discuss the gapped states, which have 
spontaneous layer polarizations with spin or valley dependent 
signs, and then ungapped states with lowered rotational sym- 
metry. 

In a continuum model, the energy of the gapped states is 
minimized when half of the spin- valley components are po- 
larized toward one layer and half to the other. 8 ^ We consider 
only states of this type, which are favored over other closely 
related states by electrostatic interactions. In a lattice model 
there is a clear distinction and an energy difference between 
states with opposite layer polarizations for different valleys, 
which have either an anomalous Hall (AH) effect or a spin 
Hall (SH) effect, and states with opposite layer polarization 
for opposite spins (LAF), which form an antiferromagnetic 
state. The three types of ABC trilayer gapped states that have 
no overall layer polarization are compared in Table I. In our 
mean-field calculations anomalous Hall and spin Hall states 
have the same energy. 

We define the condensation energy of the LAF and AH/SH 
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TABLE I: Upper panel: Mean-field theory properties of the three 
balanced-charge-density gapped states. Each of these states has two 
of the four valley/spin flavors polarized towards the top layer and 
two toward the bottom layer. Each polarized flavor contributes three 
quantized e 2 /h units to the Hall conductivity with a sign that de- 
pends on both valley and layer polarization; the continuum model as- 
signments can be retained in a lattice model because the momentum 
space Berry curvatures are strongly localized near Brillouin-zone 
corners. Middle panel: The density transfer from one outer layer 
to the other Ani for each valley- spin degree of freedom in units of 
10 11 cm~ 2 . (This density scale corresponds to ^ 1.3 • 10 -5 electrons 
per carbon atom.) The total amount of charge transferred per valley- 
spin is larger in the more stable LAF configuration than the AH or 
SH configurations. A gap is the energy gap in meV. The condensation 
energies AE cond = E gapped — E gap \ ess shown are differences between 
the ordered gapped and gapless normal phases in units of 10 -7 eV 
per carbon atom. The gapless normal state energies have been ob- 
tained from a self-consistent calculation starting from the band or- 
bital seed. The anomalous Hall and spin Hall states have the same 
energy in mean-field theory. Lower panel: Differences in total en- 
ergy between LAF and AH/SH states, AE = E LAF -E AH / SH , in units 
of 10- 9 eV per carbon atom. The exchange energy difference per 
spin/valley is separated into an intravalley (AE% K ) and an intervalley 
(AE% K ) contribution. Note that the total exchange energy difference 
satisfies AEx = 4(AE% K + AE% K ). Intervalley exchange, normally 
neglected in continuum models, makes a substantial contribution to 
the energy difference between LAF and anomalous Hall states. 



gapped states as their energy relative to the ground state en- 
ergy of the unbroken symmetry states. The unbroken symme- 
try state energy is determined by carrying out self-consistent 
mean-field calculations that are seeded by the non-interacting 
electron ground state. We find that the condensation energies 
for the ordered states are ~ 10 -7 eV per carbon atom. The 
condensation energy is approximately three times smaller than 
the product of the energy gap A gap and the charge transferred 
between layers within individual spins and valleys Ani . 

The condensation energy scales are approximately five 
times larger than those obtained for bilayer graphene 8 with 
similar approximations, presumably because the crossing 
bands are even flatter in the trilayer case, increasing the role of 



interactions. The band gaps we calculate and present in Table 
I are roughly ten times larger than the the spontaneous gap val- 
ues rsj 6 meV estimated from transport measurements in ABC 
trilayer graphene 16 , and between 1.6 to 2 times larger than the 
gaps (~ 30 meV) obtained for the bilayer graphene using the 
same value of £ r = 4. [ 8 ] This could suggest that screening 
effects are underestimated by this value of £ r , or that other in- 
teraction effects than are absent in mean-field-theory play an 
essential role. 

Experimentally the ratio of trilayer to bilayer gaps is ~ 2.5, 
close to the ratio we obtain. This suggests that the choice 
£ r = 4 quantitatively overestimates exchange effects in both 
cases. The discrepancy between theory and experiment could, 
however, be due in part to the unfavorable influence of disor- 
der in experimental samples, and also in part to inaccuracies 
in our band structure model. From the results in Table I we 
can observe, for example, that the gapped states are strongly 
suppressed when 72 and 73 have opposite signs, separating the 
Fermi points of the three Dirac cones. On the other hand, 
when 72 and 73 have the same signj^the overall effect is that 
of restoring the approximate circular symmetry of the bands, 
enhancing the chances for a gapped phase. 

Our calculations find that the nematic broken symmetry 
state is not stable in ABC trilayers. When we iterate the 
Hartree-Fock equations starting from a nematic seed, lat- 
tice rotational symmetry is restored at convergence. In both 
72 = 73 = 0, and the more realistic 72, 73 7^ case, the same 
unbroken symmetry state with three band crossing points is 
reached for self-consistent calculations starting from either 
nematic or band seeds. This result is different from the one 
obtained in the graphene bilayer case, in which the same type 
of calculation yields a stable gapless state which lowers the 
crystal's rotational symmetry giving rise to a nematic orderPSl 

The gapped solution of the Hartree-Fock equations lowers 
the total energy of the system by avoiding rapid in-plane xy 
rotation of the sublattice pseudospin direction near the band 
crossing point. 7 The gapless nematic phase lowers the total 
energy of the system by reducing the wavevector dependence 
of inter- site phase differences and introduces an anisotropic 
renormalization of the band velocity. The competition be- 
tween the two broken symmetry phases depends on how much 
energy can be gained by reshaping the quasiparticle bands in 
two different ways. Fig. 2 shows the band structures ob- 
tained from self-consistent calculations with gapped and ne- 
matic seeds. The remote hopping terms introduce a triangular 
distortion in the shape of the bands near the Fermi energy, 
but do not greatly influence gapped state properties. These 
distortions will have important consequences for the elec- 
tronic properties of doped ABC trilayers. It is noteworthy 
that the gapless phase within the minimal model develops a 
three Fermi point structure due to electron interactions alone, 
although it does not lower rotational symmetry. The inclusion 
of the 72, 73 terms also plays a role in determining the orienta- 
tion of the triangular deformation the bands undergo near the 
Dirac points in the gapless state. 

Motivated by uncertainty in the values of the remote hop- 
ping process parameters, we have performed self-consistent 
calculations over a range of values of the 72, 73 and Jn pa- 
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FIG. 3: (Color online) Energy difference between gapped and gap- 
less states AE = Eg apped - E gapless as a function of y 2 , 73 and y N . 
For Yn = the lattice symmetry of the gapless state is not lowered 
by interactions. The vertical black solid lines indicate the hopping 
parameters 72 = 0.01 eV and 73 = 0.3 eV that best approximate the 
band structure predicted by ab initio DFT calculations. For strong 
remote hopping processes, the gap in the gapped state closes pro- 
gressively and the energy difference between gapped and ungapped 
states is reduced progressively. 



rameters. The dependence of the energy difference between 
the interaction driven gapped and gapless states on model pa- 
rameters is summarized in Fig. 3. We find that the gapped 
phase almost always has a lower total energy than the gapless 
phase. However, as expected, when the remote hopping pro- 
cesses are stronger, the difference in the total energy between 
the gapped and gapless phases become smaller. Fig. 3 shows 
that the occurrence of the gapped phase relies on the principal 
intra and interlayer processes whose strength is defined by 70 
and 71 , and by the flatness of the crossing between conduction 
and valence bands that their dominance implies. 



IV. DISCUSSION 

We have used Hartree-Fock mean-field-theory calculation 
to demonstrate that electron interactions can lead to ordered 
phases in ABC trilayer graphene, provided that the strengths 
of remote hopping process in this two-dimensional crystal are 
close to current estimates. In ABC trilayers, bands near the 
Dirac point are strongly influenced by the 72 parameter over 
energy scales of ~ 20 meV, compared to the ~ 1 meV scale 
over which analogous processes play a role in AB bilayers. 
The physics of their interplay with interactions is therefore 
less likely to be distorted by disorder. Remote hopping pro- 
cesses in ABC trilayers can be important in fixing the shape 
of the energy bands near the Fermi level. We have shown 
that the gapped broken symmetry phases are nevertheless pre- 
ferred energetically over gapless states for a wide range of re- 
mote hopping parameters. We find that our gapless solutions 



do not lower the crystal symmetry, although they do generally 
lead to the formation of a triple Dirac point at the vertices of an 
equilateral triangle. The nematic phase that would break the 
triangular crystal symmetry is not stable. When remote hop- 
ping processes are included the location of the Dirac points is 
fixed by the 72 process. 

There are three distinct gapped states which have very sim- 
ilar energies. Among these the anomalous Hall and spin Hall 
(AH and SH) states have the same energy within mean-field 
theory, whereas the layer antiferromagnet state is distinct and 
is favored by inter- valley exchange because electrons with the 
same spin state have the same sense of layer polarization. 
We find that the difference in total energy between LAF and 
AH/SH is two orders of magnitude smaller than the conden- 
sation energy of either state. These states should therefore, in 
our view, be considered as close cousins. It seems likely that 
real samples are likely to be found in configurations in which 
several phases are present separated by domain walls. An ex- 
ternal magnetic field which favors anomalous Hall states, at 
least at finite carrier densities, likely can be used to manipu- 
late the domain structure. 

Using the Hartree-Fock approximation and reducing inter- 
action strengths by a factor of £ r = 4, we find that ABC tri- 
layer graphene has a substantial interaction driven gap of the 
order of 65 meV. The size of the gap is sensitive to the choice 
we have made for the £ r parameter, which we have chosen to 
mimic exchange interaction renormalization parameters that 
are used in ab initio hybrid-density-functional calculations. 
Band structure effects can reduce the size of the gap, sub- 
stantially so when 72 is assigned a negative value. For fa- 
vorable parameters the gaps we find are approximately twice 
as large as than those predicted values for bilayer graphene 
using corresponding approximations. The theoretical gaps are 
therefore very much larger than initial estimates of a sponta- 
neous band gap from ABC trilayer experiments which suggest 
a value ~ 6 me\4^ The discrepancy is certainly due in part 
to disorder and inhomogeneity which reduces the gaps of ex- 
perimental systems below ideal values, but could also reflect 
a theoretical overestimate. We note in this connection that 
ABC trilayer graphene samples generally have poorer quality 
than bilayers. This difference could be due to lower effective- 
ness of the current annealing procedure routinely applied to 
suspended graphene single or multilayer samples. Future ex- 
perimental work may establish a higher lower bound for the 
trilayer graphene gap. 
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